Non-Markovian dynamics without using quantum trajectory 
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Open quantum system interacting with structured environment is important and manifests non- 
Markovian behavior, which was conventionally studied using quantum trajectory stochastic method. 
In this paper, by dividing the effects of the environment into two parts, we propose a determin- 
istic method without using quantum trajectory. This method is more efficient and accurate than 
stochastic method in most Markovian and non-Markovian cases. We also extend this method to the 
generalized Lindblad master equation. 

PACS numbers: 03.65.Yz, 42.50.Lc 



When an open quantum system interacts with environ- 
ment, it experiences decoherence and dissipation which 
lead to loss of information. Such open quantum system is 
depicted by a reduced density matrix which shows non- 
unitary evolution. On the other hand, the environment is 
classified as Markovian with no memory effect, and non- 
Markovian with memory effect. In Markovian case, since 
there is no memory effect, the quantum trajectory based 
Monte Carlo wave function (MCWF) method [1-3] and 
quantum state diffusion (QSD) method [4, 5] are applied. 
However, in non-Markovian case, due to memory effect, 
the information lost by the system during the interaction 
with the environment will come back to the system in a 
later time and so shows much more complicated behav- 
iors than Markovian case. 

Non-Markovian systems are important for their ap- 
plications to many fields of physics, such as quan- 
tum information processing [6, 7], quantum optics [8], 
solid state physics [9], and chemical physics [10]. Re- 
cently, non-Markovian behaviors have also been studied 
in biomolecules where the molecules are embedded in 
a solvent and/or in a protein environment [11]. Since 
there is no true pure state quantum trajectory due to 
the memory effect [12], the quantum trajectory based 
Markovian methods do not work. Thus, doubled Hilbert 
space (DHS) method [13], triple Hilbert space (THS) 
method [14], non-Markovian QSD method [15, 16], and 
non-Markovian quantum jump (NMQJ) method [17, 18] 
are proposed to solve the non-Markovian dynamics of the 
system where the memory effect is taken into account. 
However, in order to obtain high accuracy, all these meth- 
ods, which are based on stochastic simulations, need to 
fulfill a large number of realizations and is very time- 
consuming. So, new methods which are more efficient 
and accurate are highly desired. 

In this paper, a deterministic method without us- 
ing quantum trajectory is proposed to solve the non- 
Markovian dynamics. The influence of the environment 
on the system is divided into two parts, i.e., the non- 



unitary evolution of the states and the probability flow 
between these states. Moreover, we also extend this 
approach to the generalized Lindblad master equation 
which can deal with some strong coupling cases [19]. The 
algorithm and numerical efficiency are given, which show 
that our method is more efficient and accurate than those 
based on stochastic simulation in most Markovian and 
non-Markovian cases. 

The dynamics of the non-Markovian system is gov- 
erned by the following master equation [8] 



p(t) =^\H s ,p{t)]+Y,lAt)C 3 {t)p{t)C){t) 
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where H s is the system Hamiltonian including the Lamb 
shift, Cj(t) are the jump operators which induce changes 
[e.g., jump from state ip a (t) to ip a '(t) i.e., \tp a >(t)) = 
Cj(t)\4> a {t))/ \\Cj{t)\ip a (t))\\ ] in the system, and Jj(t) 
are the decay rates which may take negative values for 
some time intervals. The reduced density matrix can be 
written as [17] 



P (t) = J2p<*{t)\Mt))(Mt)\ 
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where p a (t) is the probability of the system being in the 
state \ip a (i)) at time t. Further, it should be pointed 
out that the effective number of the states N e ff is de- 
termined by Cj(t)'s [18], J2a=i P<*(t) = 1 and that the 
state \i/} a (t)) is normalized. 

To solve the dynamics of the system, one should know 
the time evolution of \tp a (t)) and its probability p a (t). In 
our method, the time evolution of the state \ip a (t)) is the 
same as that in NMQJ [17]. In NMQJ, the probability 
p a (t) is calculated in a stochastic way by using quan- 
tum trajectory to N ensemble members. In our method, 
however, the evolution of probability p a (t) is given in a 
deterministic way: 

Mt) = -J2 r « (*) + E ' T «> (*)*<»' (*)> ( 3 ) 
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where T J a (t) — jj(t)\\Cj(t) \ip a (t))\\ and ^'represents 

,j) 

the summation over all the pairs (a',j) satisfying 
\i/j a (t)) = Cj{t) \ijj a ,{t))/\\C 3 {t) \ip a >{t))\\. One finds that 
the probability of the state, p a (t), changes via the mech- 
anism of jumps for "out" (a — > a') and "in" (a' — > a), 
respectively. 

The numerical simulation corresponding to Eq. (3) is 
straightforward: 

Pa {t+St) = Pa (t)-5tJ2 T Ut)Pc*(t)+5t E H'WM*). 

3 (a',j) 

(4) 

Note that there is no stochastic noise and no need to 
consider the sign of the decay rate during the simula- 
tion. Additionally, the p Q (i)'s in our method do repre- 
sent the probability of the system actually being in the 
corresponding pure state ensemble. 

Consider a particular transition: \tp a i(t)) = 
Cj(t)\ifj a (t))/ \\Cj(t)\ip a (t)}\\, then the corresponding 
probability change takes the form: 



p a (t + 8t) =p a (t)-St Pa (t)Ti(t), 
pctt + St) = Pa ,{t) + 6tp a (t)T j a (t). 



(5) 



When the decay rate jj(t) is positive or negative, the 
probability flow is from \ip a (t)) to \ip a '(t)) or reversed. 
This has been mentioned in Rcf.[17]. However, it is more 
explicit in our method. From Eq. (5), it is clear that, in 
the negative decay region, the amount of probability flow 
only depends on the target state and the probability of 
the system being in the target state. This is similar to 
the situation in NMQJ [17], where the jump probability 
in the negative decay region is proportional to the num- 
ber of particles in the target state. These indicate that 
the trajectory of a particle in NMQJ can not be inter- 
preted as true trajectory since the jump process depends 
on the status of other particles in the system. Because 
true pure state quantum trajectories do not exist in the 
non-Markovian dynamics [12], it is not necessary to cal- 
culate p a (i) in a stochastic way. 

Next, we extend our method to the recently proposed 
generalized Lindblad master equation which can solve the 
dynamics of some highly non-Markovian systems[19], 



d 
dt 
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where i,j = 1, 2, • • ■ , n, Hi are any Hcrmitian operators 
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and RV are any system operators. It should be indicated 



that p(t) = E Pitt)- 
»=i 

The ith density matrix is decomposed as: 



E 
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where N % e ^ is determined in the same way in Eq. (2) 
by taking all the jump operators R z J's and all the states 
^(tys in each Pj(t) into consideration. 

The evolution of state \ipf{t)) is governed by the non- 
linear differential equation [21] 



(8) 



where G(ipf)(t) 
fE||^lV>f(*)}f 
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By combining Eqs. (6), (7), 



(8) and noting that \^}(t)) (^j(t)\, \tf(t)) <V*C*)|i ■ "> 
i£" ff (t) are linearly independent, the 

evolution of pf (t) is given by 

#(t) = -£lM(t)+ E ' T 2a>P?(t)> ( 9 ) 
3" (p,",a') 

where T l J a = | V^* (*) ) 1 1 2 an <4 E ' represents the 
summation over all the pairs (J, v, a') satisfying \ipf(t)) = 
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It can be easily seen that 



by setting n = 1 and taking the decay rates j(t) into the 
equation, Eq. (9) degenerates to Eq. (3). 

Example 1: Detuned Jaynes-Cummings model- 
Consider a system with a two-level atom in a detuned 
damped cavity, which is governed by the time convolu- 
tionless master equation [8] 

p(t)=-^S(t){<r+<r-,P(t)} 

(10) 

The spectral density of the cavity is supposed to be of 

A 2 

Lorentzian profile, i.e., J(u) = 27r[(aJo _A_^) 2 +A 2 ] > where 
A = ujq — uj c is the detuning between the cavity mode 
and the atom. To second order approximation, the 
Lamb shift and the decay rate take the form [8] S(t) = 



70 A A 
A 2 +A 2 
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{l-c- A *[cos(A£) + ^sin(A£)]}, 7 (f) 
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e _A *[cos(At) — -j- sin(Ai)]}. In this model, there is only 
one jump operator C = er_ = \g) (e|, which is a lower- 
ing operator. We assume that p(0) = |^i(0)) (^i(0)| and 
choose IV'i(O)) = (4 |e)+3 \g))/5. Acting the jump opera- 
tor on the state |?/>i(0)), we get \ip 2 {0)) — \g)- According 
to Eq. (4), at time t + St, the probabilities become 



Pi(t 

P2(t 



St) 
St) 



Pl(t) 

P2(t) 



St Pl (t)T\(t), 
St Pl (t)Tl(t), 



(11) 



where T\(t) = 7 (t) \(e \ Vi(*)>| 2 ■ 

In this example, p ee (t) is proportional to the energy of 
the system and p 2 (t) represents the probability for one 
photon being in the environment. Although pi(t) and 
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P2(t) can be solved analytically, in order to illustrate our 
method, we use Eq. (11) to do the simulation. The pa- 
rameters are chosen as A = 12A,7oA = 4, XSt = 0.005. 
Figure 1 (a) shows explicitly the reversal of the proba- 
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FIG. 1: (color online) Dynamics of detuned Jaynes- 
Cummings model. The initial state is \ipi(0)) = (4 \e) + 
3 \g})/5 and the parameters are A = 12A, 7A = 4, XSt = 0.005. 
(a) The probabilities for the system in states \tjJi{t)) an d 
\ip2(t)). (b) The population of the excited state p ee (ini- 
tially higher line) and the absolute value of the coherence p eg 
(initially lower line) with three methods: analytic (red solid 
curve), our method (blue long-dashed curve) and NMQJ (with 
N = 10 4 particles in the system, green dash-dot curve). 

bility flow. We can see from Fig. 1 (a) and (b) that 
when the probability flow gets reversed, the energy and 
coherence of the atom increase. These show explicitly the 
memory effect that the reduced system restores the in- 
formation lost earlier. In Fig. 1 (b), the result of NMQJ 
(with N = 10 4 particles in the system) is also given, 
which shows that our method is more accurate. 

Example 2: Application to generalized Lindblad mas- 
ter equation - To illustrate our method for this kind of 
equation, we consider a two-state system coupled to an 
environment consisting of two energy bands, each with 
a finite number of evenly spaced levels. This may be 
viewed as a spin coupled to a single molecule or a single 
particle quantum dot [20] . By using time-convolutionless 
projection operator technique, to the second order, the 
generalized Lindblad master equation takes the form [21] 

d f l 

—P! = / dtlfl(t - tl)[2^/l(J + p 2 (T~ - 7 2 {<7 + (7~, / 9l}], 

at J 

—P2 = / dhh(t - ii)[272er~picr + - 7i{er~er + , p 2 }], 
at J 

(12) 

where Jih(t — ti),(i = 1,2), is the environment corre- 
lation function with h(t) = SC 2^(Set%)^ > wnere Se is the 
width of the upper and lower energy bands. The reduced 
density matrix for the system is given by p = p\ + p2- 
We assume that pi(0) = |e) (e| and p2(0) = 0. The pa- 
rameters are chosen as Se = 0.31 and 71 = 72 = 1. In 
Fig. 2 we compare the results of our method, analytical 
solution and Monte Carlo simulation which is based on 
the unraveling of the master equation (with N = 10 tra- 
jectories) [21]. Apparently, our method is more accurate 
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FIG. 2: (color online) A two-state system coupled to an envi- 
ronment consisting of two energy bands. Comparison of our 
method (blue long-dashed curve) and Monte Carlo simulation 
(with N = 10 4 trajectories, green dash-dot curve) to analyt- 
ical result (red solid curve). The parameters are Se = 0.31, 
71 = 72 = 1 and time step St = 0.01. 



than Monte Carlo simulation method. 

According to Eqs. (4), we only need to calculate N e ff 
states and change the probabilities dctcrministically. The 
time cost is almost determined by the calculation of N e f f 
states. However, the evolution of N e ff states is indepen- 
dent with each other, so we can calculate them parallelly. 
In addition, if the jump operators can be represented by 
sparse matrixes, we only need to calculate the evolution 
of the states appearing in the decomposition of p(0) and 
use the jump operators to obtain other states. Moreover, 
since the sign of the decay rate makes no difference dur- 
ing the simulation, in non-Markovian case, our method 
is as efficient as it behaves in Markovian case. 

Similar to our method, the NMQJ method [17, 18] 
needs to calculate N e ff states. However, in addition to 
that, NMQJ has to consider the sign of the decay rates 
and generate N random numbers (TV 3> N e ff) to decide 
the jump process at each time step St. Apparently, our 
method is more efficient than NMQJ in any case. 

In Markovian case, the MCWF [1] and QSD [4] method 
need to realize a large number of trajectories for every 
state appearing in the decomposition of p(0). When 
the number of these trajectories is larger than N e ff, 
which is always the case, our method is more efficient 
than them. In non-Markovian case, the DHS method 
[13], THS method [14] and non-Markovian QSD method 
[15, 16] all introduce additional cost for computational 
efficiency compared to MCWF or QSD. However, in non- 
Markovian case, our method is as efficient as it behaves in 
Markovian case. Thus, when the number of these trajec- 
tories is larger than N e f f , our method is obviously more 
efficient than them, too. 

As for the accuracy, since there is no statistical noise 
in our method and the error caused by finite time step 
St is the same, compared with all the methods based 
on stochastic simulation, our method is more accurate. 
Actually, our method is the limit case when the number 
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of realizations in the stochastic based methods tends to 
infinite. 

In conclusion, by dividing the influence of the environ- 
ment on the system into two parts, i.e., the non- unitary 
evolution of these states and the probability flow between 
them, we propose a deterministic method to solve the 
non-Makovian dynamics. Compared with the method 
based on stochastic simulation, our method has advan- 



tages in efficiency and accuracy. Additionally, we ex- 
tended this approach to the generalized Lindblad master 
equation , which is useful to solve the dynamics of some 
highly non-Markovian systems. 
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